L9: Rasters and Raster Algreba

Reading Rasters, Raster Attributes, Raster Algebra

Bogdan G. Popescu
bogdan.popescu@johncabot.edu

John Cabot University

Today

Todays agenda includes:

  • becoming familiar with others types of data in R:

  • raster layers

  • examining spatial and non-spatial properties of vectors

We will use sf, stars, and units packages.

Reminder: Spatial Data Classes in R

  • Spatial entries in GIS can be represented as a vector or as a raster

  • In the previous lecture, we examine vectors and geometric operations with vectors.

Raster

Rasters are matrices or an array representing a rectangular area on the surface of the earth.

They contain elements such as :

  • Origin (\(x_{max}\), \(y_{min}\)) and resolution (\(\delta_{x}\), \(\delta_{y}\))
  • Coordinate Reference System
  • Values
  • Dimensions (rows, column, layers)
  • Extent: \(x_{min}\), \(y_{min}\), \(x_{max}\), \(y_{max}\)

Raster

Raster with the Same Extent but Different Resolutions

Raster file formats

  • “Simple” rasters
    • GeoTIFF (.tif)
    • Erdas Imagine Image (.img)
  • “Complex” rasters (>3D and / or metadata)
    • HDF (.hdf, .he5)
    • NetCDF (.nc)

Using R packages for rasters

There is one important package that we will work with: stars

terra is also another important package developed by the same authors

Both contain classes, and extensive functions for working with rasters in R

We will mostly work with stars

The stars package

Unlike terra and raster, stars is very well integrated with sf - the package that we have used for vector analysis.

Example of Raster: Luminosity

  • A typical example of raster is satelite luminosity

  • Within economics, researchers use satelite luminosity as a proxy for economic activity

Luminosity in 2013

Working with Luminosity

You can download satelite luminosity by going to https://www.ncei.noaa.gov/products/dmsp-operational-linescan-system


Working with Luminosity

The next step is to scroll down to products


Working with Luminosity

The next step is to scroll down to products


Working with Luminosity

You can download these files by clicking on the links


Working with Luminosity

Note that these are large files: Every zip file will range from 355 MB to 411.1 MB

The files are cloud-free composites

In cases where two satellites were collecting data - two composites were produced.

Working with Luminosity


Working with Luminosity

Note that these are large files: Every zip file will range from 355 MB to 411.1 MB

The files are cloud-free composites

In cases where two satellites were collecting data - two composites were produced.

The products are 30 arc second grids, spanning -180 to 180 degrees longitude and -65 to 75 degrees latitude.

Each composite set is named with the satellite and the year (F121995 is from DMSP satellite number F12 for the year 1995)

Working with Luminosity

Three types of images are available:

  • F1?YYYY_v4b_cf_cvg.tif: Cloud-free coverages tally the total number of observations that went into each 30 arc second grid cell
  • F1?YYYY_v4b_avg_vis.tif Raw avg_vis contains the average of the visible band digital number values with no further filtering.
  • F1?YYYY_v4b_stable_lights.avg_vis.tif The cleaned up avg_vis contains the lights from cities, towns, and other sites with persistent lighting, including gas flares.

Working with Luminosity

The DMSP/OLS resolution is approximately 1000m

It comprises six different DMSP satellites:

  • F10 (1992–1994)
  • F12 (1994–1999)
  • F14 (1997–2003)
  • F15 (2000–2007)
  • F16 (2004–2009)
  • F18 (2010–2013)

Working with Luminosity

After downloading, your folder should look like below:


Working with Luminosity

You download all the relevant files for this exercise below:

Working with Luminosity

These are the steps to have the data ready for plotting in ggplot

#Step1: Read country borders
borders <- st_read(dsn="./data/boundaries/all-country-boundaries.shp", quiet = TRUE)
#Step2: Creating a folder called
target_dir <- "./data/luminosity/F101992"
#Step3: Create the target directory if it doesn't exist
if (!dir.exists(target_dir)) {
  dir.create(target_dir)
}
#Step4: Unzip the file
untar("./data/luminosity/F101992.v4.tar", exdir = target_dir)
#Step5: Unzip the file further
R.utils::gunzip("./data/luminosity/F101992/F101992.v4b_web.avg_vis.tif.gz", remove = FALSE, overwrite=TRUE)
#Step6: Read the raster
f1992<-read_stars("./data/luminosity/F101992/F101992.v4b_web.avg_vis.tif")
#Step7: Make sure the two files have the same CRS
st_crs(f1992)<-st_crs(borders)
#Step8: Reduce the resolution of the raster
f1992b<-st_warp(f1992, st_as_stars(st_bbox(f1992), dx = 0.7))
#Step9: Rename the raster
names(f1992b)<-"lights_1992"

fig<-ggplot()+
  geom_stars(data=f1992b)+
  scale_fill_gradientn(colours=c("black","white"))+
  #scale_fill_viridis_c(option = "magma",begin = 0.1)+
  geom_sf(data=borders, linewidth = 0.1, fill = NA, color = "white", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

fig

Luminosity in 1992

Luminosity in 1992 in Europe

Raster Values and Properties

The print method gives a summary of the raster’s properties.

print(f1992b)
stars object with 2 dimensions and 1 attribute
attribute(s):
             Min. 1st Qu. Median     Mean 3rd Qu. Max. NA's
lights_1992     0       2      3 3.384738       4   63  515
dimension(s):
  from  to offset delta refsys x/y
x    1 515   -180   0.7 WGS 84 [x]
y    1 201     75  -0.7 WGS 84 [y]

The class function allows us to see the type of object that we’re dealing with

class(f1992b)
[1] "stars"

Raster Values and Properties

Fundamentally, a stars object is a collection of matrix or array objects.

This is along with additional properties of the dimensions, such as dimension names, Coordinate Reference Systems (CRS).

We can display the structure of a stars object using str.

str(f1992b)
List of 1
 $ lights_1992: num [1:515, 1:201] 0 0 0 0 0 0 0 0 0 0 ...
 - attr(*, "dimensions")=List of 2
  ..$ x:List of 7
  .. ..$ from  : num 1
  .. ..$ to    : num 515
  .. ..$ offset: num -180
  .. ..$ delta : num 0.7
  .. ..$ refsys:List of 2
  .. .. ..$ input: chr "WGS 84"
  .. .. ..$ wkt  : chr "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223"| __truncated__
  .. .. ..- attr(*, "class")= chr "crs"
  .. ..$ point : logi NA
  .. ..$ values: NULL
  .. ..- attr(*, "class")= chr "dimension"
  ..$ y:List of 7
  .. ..$ from  : num 1
  .. ..$ to    : num 201
  .. ..$ offset: num 75
  .. ..$ delta : num -0.7
  .. ..$ refsys:List of 2
  .. .. ..$ input: chr "WGS 84"
  .. .. ..$ wkt  : chr "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223"| __truncated__
  .. .. ..- attr(*, "class")= chr "crs"
  .. ..$ point : logi NA
  .. ..$ values: NULL
  .. ..- attr(*, "class")= chr "dimension"
  ..- attr(*, "raster")=List of 4
  .. ..$ affine     : num [1:2] 0 0
  .. ..$ dimensions : chr [1:2] "x" "y"
  .. ..$ curvilinear: logi FALSE
  .. ..$ blocksizes : NULL
  .. ..- attr(*, "class")= chr "stars_raster"
  ..- attr(*, "class")= chr "dimensions"
 - attr(*, "class")= chr "stars"

Raster Values and Properties

From this, we see that we have the f1992b object of length 1.

The 2nd line in the str printout specifies the name and contents of the first (and only) item in the list, namely its type, dimensions, and a sample of the first few values.

In our case, the first item is lights_1992, and it is a matrix with 515 rows and 201 columns

Raster Attributes and Values

To access some of the attributes of our objects, we can type:

names(f1992b)
[1] "lights_1992"

We can change names easily in the following wat:

names(f1992b)<-"lum_1992"
names(f1992b)
[1] "lum_1992"

Accessing Attributes by Name

We can access our raster’s attributes by using the list access methods

#f1992b$lights_1992
f1992b[["lum_1992"]][1:10, 1:10]
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0    0    0    0    0    0    0   10    7     9
 [2,]    0    0    0    0    0    0    0    6    7     6
 [3,]    0    0    0    0    0    0    0    5    6     9
 [4,]    0    0    0    0    0    0    0    9    6     9
 [5,]    0    0    0    0    0    0    0    8    8     8
 [6,]    0    0    0    0    0    0    0    9    9     7
 [7,]    0    0    0    0    0    0    0    8    7     8
 [8,]    0    0    0    0    0    0    0    4   10     9
 [9,]    0    0    0    0    0    0    0    6    7     9
[10,]    0    0    0    0    0    0    0    0    8     9

This is how we demonstrate that our file is a matrix

#Using names
class(f1992b$lum_1992)
[1] "matrix" "array" 
#Using index
class(f1992b[[1]])
[1] "matrix" "array"

Accessing Attributes by Name

And that it has 512 rows and 201 columns.

#Using names
dim(f1992b$lum_1992)
  x   y 
515 201 
#Using index
dim(f1992b[[1]])
  x   y 
515 201

Subset operators in R

Syntax Objects Returns
x[i] vector, list Subset i
x[i, j] data.frame, matrix Subset i, j
x[i, j, k] array Subset i, j, k
x[[i]] vectors, lists Single element i
x$i data.frame, list Single column or element i
x@n S4 objects Slot n

Dimensions and spatial properties

We can use the following function to get at important raster properties

# Number of rows
nrow(f1992b)
  x 
515 
# Number of columns
ncol(f1992b)
  y 
201 
# All dimensions
dim(f1992b)
  x   y 
515 201 
# Extent of a raster
st_bbox(f1992b)
      xmin       ymin       xmax       ymax 
-180.00417  -65.69583  180.49583   75.00417

Dimensions and spatial properties

The other properties of rasters can be extracted in the following way:

st_dimensions(f1992b)$x$offset  ## x-axis origin
[1] -180.0042
st_dimensions(f1992b)$y$offset  ## x-axis origin
[1] 75.00417
st_dimensions(f1992b)$x$delta  ## x-axis resolution
[1] 0.7
st_dimensions(f1992b)$y$delta  ## y-axis resolution
[1] -0.7

Note that the resolution is separate for x and y.

If the (absolute) resolution is euqal for both, raster pixels are square.

When the x and y resolutions are unequal, raster pixels are non-square rectangles.

Dimensions and spatial properties

We can find out the crs of our raster using the st_crs function

st_crs(f1992b)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

Accessing Raster Values

As indicated, raster values can be extracted directly, as a matrix or as an array.

For example:

#Printing only the first 10 rows and first 10 columns
f1992b$lum_1992[1:10, 1:10]
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0    0    0    0    0    0    0   10    7     9
 [2,]    0    0    0    0    0    0    0    6    7     6
 [3,]    0    0    0    0    0    0    0    5    6     9
 [4,]    0    0    0    0    0    0    0    9    6     9
 [5,]    0    0    0    0    0    0    0    8    8     8
 [6,]    0    0    0    0    0    0    0    9    9     7
 [7,]    0    0    0    0    0    0    0    8    7     8
 [8,]    0    0    0    0    0    0    0    4   10     9
 [9,]    0    0    0    0    0    0    0    6    7     9
[10,]    0    0    0    0    0    0    0    0    8     9

Accessing Raster Values

A histogram will help us identify the raster values distribution

#Turning the object into a dataframe
df<-st_as_sf(f1992b[1], as_points = FALSE, merge = FALSE)

#Creating a histogram
ggplot(df, aes(x = lum_1992)) +
  geom_histogram(binwidth = 10, color = "white") +
  ggtitle("Histogram of lum_1992") +
  xlab("lum_1992") +
  ylab("Frequency")+
  theme_bw()

Accessing Raster Values

We thus see that most of values around the world are 0.

This is because 71 percent of the Earth’s surface is water-covered.

The important part is that by transforming the stars object into a dataframe (grid), we can perform the usual operations dedicated to dataframes.

For example: glimpse:

library(dplyr)
glimpse(df)
Rows: 103,000
Columns: 2
$ lum_1992 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ geometry <POLYGON [°]> POLYGON ((-180.0042 75.0041..., POLYGON ((-179.3042 7…

Accessing Raster Values

We can also calculate mean, sd, min, and max

mean(df$lum_1992, na.rm=T)
[1] 3.384738
sd(df$lum_1992, na.rm=T)
[1] 3.161684
min(df$lum_1992, na.rm=T)
[1] 0
max(df$lum_1992, na.rm=T)
[1] 63

Mapping the World

Mapping Lights in Central Rome

Let’s focus on Italy.

You can download the Rome central area polygon

You can also download the main streets for the central area

Make sure you place them in the relevant folder.

Mapping Lights in Central Rome

Let’s focus on Italy

#Step1: Reading the central polygon
relev_poly<-st_read(dsn="./data/rome_center/central_polygon.shp", quiet = T)
#Step2: Reading roads
gis_osm_roads <- st_read(dsn="./data/selection.gdb", layer="gis_osm_roads", quiet = T)
#Step3: Extracting the bix from the central polygon
ext2<-st_bbox(relev_poly)
#Step4: Ensuring the same crs for the polygon as for the raster
st_crs(relev_poly)<-st_crs(f1992)
#Step5: Cropping the raster
rc <- st_crop(f1992, relev_poly)

#Step6: Mapping everything
ggplot()+
  geom_stars(data = rc)+
  geom_sf(data = relev_poly, color="white", fill=NA)+
  geom_sf(data = gis_osm_roads, color="white", linewidth=0.05)+
  theme_bw()+
  scale_fill_gradientn(colours=c("black","white"))

Mapping Lights in Central Rome

Let’s focus on Italy

Mapping Lights in Central Rome

Every grid has a value:

#Turning the object into a dataframe
rc2<-st_as_sf(rc[1], as_points = FALSE, merge = FALSE)
names(rc2)[1]<-"lum_1992"
#Step6: Mapping everything
ggplot()+
  geom_sf(data = rc2, aes(fill=lum_1992), linewidth=0)+
  geom_sf(data = relev_poly, color="white", fill=NA)+
  geom_sf(data = gis_osm_roads, color="grey50", linewidth=0.05)+
  theme_bw()+
  scale_fill_gradientn(colours=c("black","white"))+
  geom_sf_text(data = rc2, 
  aes(label = lum_1992, geometry = geometry),
  color = ifelse(rc2$lum_1992 < 61, "white", "black"),
  #aes(label = lum_1992, geometry = geometry, color = ifelse(lum_1992 < 61, "white", "black")),
  size = 1.5)

Mapping Lights in Central Rome

Every grid has a value:

Mapping Lights in Central Rome

Let’s zoom in:

#Turning the object into a dataframe
rc2<-st_as_sf(rc[1], as_points = FALSE, merge = FALSE)
names(rc2)[1]<-"lum_1992"
#Step6: Mapping everything
ggplot()+
  geom_sf(data = rc2, aes(fill=lum_1992), linewidth=0)+
  geom_sf(data = relev_poly, color="white", fill=NA)+
  geom_sf(data = gis_osm_roads, color="grey50", linewidth=0.05)+
  theme_bw()+
  scale_fill_gradientn(colours=c("black","white"))+
  geom_sf_text(data = rc2, 
  aes(label = lum_1992, geometry = geometry),
  color = ifelse(rc2$lum_1992 < 61, "white", "black"),
  #aes(label = lum_1992, geometry = geometry, color = ifelse(lum_1992 < 61, "white", "black")),
  size = 4)+
      coord_sf(xlim = c(12.45, 12.60), 
           ylim = c(41.87, 41.79))

Mapping Lights in Central Rome

Let’s zoom in:

Mapping Lights in Central Rome

As you might imagine, we can use different colors for these values:

Mapping Lights in Central Rome

As you might imagine, we can use different colors for these values:

Interractive Maps with mapview

library(mapview)
mapview(rc2)

Replacing Values

We can easily replace the values inside the raster (grid) by using simple dataframe operations.

For example, let’s replace all the values that have 63 with NA:

rc3<-rc2
rc3$lum_1992[rc3$lum_1992==63]<-NA
pic<-ggplot()+
  geom_sf(data = rc3, aes(fill=lum_1992), linewidth=0)+
  geom_sf(data = relev_poly, color="white", fill=NA)+
  geom_sf(data = gis_osm_roads, color="grey50", linewidth=0.05)+
  theme_bw()+
  scale_fill_viridis_c(option = "magma",begin = 0.1)+
  #scale_fill_gradientn(colours=c("black","white"))+
  geom_sf_text(data = rc2, 
  aes(label = lum_1992, geometry = geometry),
  color = ifelse(rc3$lum_1992 < 61, "white", "black"),
  #aes(label = lum_1992, geometry = geometry, color = ifelse(lum_1992 < 61, "white", "black")),
  size = 4)+
      coord_sf(xlim = c(12.45, 12.60), 
           ylim = c(41.87, 41.79))
pic

Replacing Values

We can easily replace the values inside the raster (grid) by using simple dataframe operations.

For example, let’s replace all the values that have 63 with NA:

Replacing Values

We can easily replace the values inside the raster (grid) by using simple dataframe operations.

For example, let’s replace all the values that have 63 with NA:

Let us now make the NA grids transparent

rc3<-rc2
rc3$lum_1992[rc3$lum_1992==63]<-NA
pic<-ggplot()+
  geom_sf(data = rc3, aes(fill=lum_1992), linewidth=0)+
  geom_sf(data = relev_poly, color="white", fill=NA)+
  geom_sf(data = gis_osm_roads, color="grey50", linewidth=0.05)+
  theme_bw()+
  scale_fill_viridis_c(option = "magma",begin = 0.1,
                       na.value = "white")+
  #scale_fill_gradientn(colours=c("black","white"))+
  geom_sf_text(data = rc2, 
  aes(label = lum_1992, geometry = geometry),
  color = ifelse(rc3$lum_1992 < 61, "white", "black"),
  #aes(label = lum_1992, geometry = geometry, color = ifelse(lum_1992 < 61, "white", "black")),
  size = 4)+
      coord_sf(xlim = c(12.45, 12.60), 
           ylim = c(41.87, 41.79))
pic

Replacing Values

We can easily replace the values inside the raster (grid) by using simple dataframe operations.

For example, let’s replace all the values that have 63 with NA:

Let us now make the NA grids transparent

Writing raster to file

We can of course write the grid as a shapefile for we can turn it back to a raster.

rc3_raster = st_rasterize(rc3)
class(rc3_raster)
[1] "stars"
class(rc3)
[1] "sf"         "data.frame"
#Step1: Creating a folder called
target_dir <- "./data/output"
#Step2: Create the target directory if it doesn't exist
if (!dir.exists(target_dir)) {
  dir.create(target_dir)
}
#Step3: Writing the raster
write_stars(rc3_raster, './data/output/rc3_raster.tif', update = TRUE)
#Step4: Writing the shape files
st_write(rc3, "./data/output/rc3_grid.shp", delete_dsn = TRUE)
Deleting source `./data/output/rc3_grid.shp' using driver `ESRI Shapefile'
Writing layer `rc3_grid' to data source 
  `./data/output/rc3_grid.shp' using driver `ESRI Shapefile'
Writing 535 features with 1 fields and geometry type Polygon.

Working with More Complex Rasters

Working with lights was simple

In many other contexts, we can have more complex rasters

Temperature data from CEDA is an example of a raster that contains time series

Temperature

Temperature in 2019

Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

You can download temperature data from by going to https://archive.ceda.ac.uk


Temperature

For a start we will simply plot a date to get a better sense of the structure of this type of raster.

library(terra)
#Step1: Read the raster
b <- rast("./data/temp/cru_ts4.07.1901.2022.tmp.dat.nc")

#Step2: Extract time information
time_values <- time(b)

#Step3: Turning strings into time
date_objects <- as.Date(time_values)

#Step4: Select only the dates from the year 1901 and 2022
dates_1901 <- unique(date_objects[format(date_objects, "%Y") == "1901"])
dates_2022 <- unique(date_objects[format(date_objects, "%Y") == "2022"])

#Step5: Find the index of the selected dates in the time_values
selected_indices_1901 <- which(time_values %in% as.Date(dates_1901))
selected_indices_2022 <- which(time_values %in% as.Date(dates_2022))

#Step6: Extract raster values for the selected dates
selected_rasters_1901 <- b[[selected_indices_1901]]
selected_rasters_2022 <- b[[selected_indices_2022]]

#Step7: Extract layer names
layer_names_1901 <- names(selected_rasters_1901)
layer_names_2022 <- names(selected_rasters_2022)

Temperature

For a start we will simply plot a date to get a better sense of the structure of this type of raster.

#Step8: Find indices of layers that start with "tmp"
#tmp stands for temperature. There are other rasters within these files.
indices_tmp_1901 <- grep("^tmp", layer_names_1901)
tmp_layers_1901 <- selected_rasters_1901[[indices_tmp_1901]]
indices_tmp_2022 <- grep("^tmp", layer_names_2022)
tmp_layers_2022 <- selected_rasters_2022[[indices_tmp_2022]]

#Step9: Changing raster names to temperature
names(tmp_layers_1901)<-dates_1901
names(tmp_layers_2022)<-dates_2022

#Step10: Turning the spatraster into raster
stars_list_1901 <- lapply(tmp_layers_1901, function(raster) st_as_stars(raster::raster(raster)))
stars_list_2022 <- lapply(tmp_layers_2022, function(raster) st_as_stars(raster::raster(raster)))

#Step11: Use lapply to apply st_apply to each stars object to calculate the mean
mean_list_1901 <- lapply(stars_list_1901, function(stars_obj) st_apply(stars_obj, c("x", "y"), mean))
mean_list_2022 <- lapply(stars_list_2022, function(stars_obj) st_apply(stars_obj, c("x", "y"), mean))

#Step12: Combine the list of results into a single stars object
mean_stars_1901 <- do.call(stars::st_as_stars, mean_list_1901)
mean_stars_2022 <- do.call(stars::st_as_stars, mean_list_2022)

#Step13: Changing the name of the rasters
names(mean_stars_1901)<-"avg_tmp_1901"
names(mean_stars_2022)<-"avg_tmp_2022"

Temperature

Temperature in 1901

fig<-ggplot()+
  geom_stars(data=mean_stars_1901)+
   scale_fill_gradient2(
    low = 'blue',
    high = 'red',
    mid = 'yellow', 
    midpoint = 5,
    na.value = NA)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "white", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Temperature

Temperature in 1901

Temperature

Temperature in 2022

fig<-ggplot()+
  geom_stars(data=mean_stars_2022)+
   scale_fill_gradient2(
    low = 'blue',
    high = 'red',
    mid = 'yellow', 
    midpoint = 5,
    na.value = NA)+
  #scale_fill_gradientn(colours=c("blue","yellow", "red"),   na.value = NA)+
  #scale_fill_viridis_c(option = "magma",begin = 0.1)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "white", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Temperature

Temperature in 2022

Working with Temperature

Now plotting was to some extent straightforward.

Now what if we wanted to extract the average temperature for every single year and make a time plot?

Working with Temperature

This is where our programming skills come in handy from earlier lectures

library(terra)

# Step 1: Read the raster
b <- rast("./data/temp/cru_ts4.07.1901.2022.tmp.dat.nc")

# Step 2: Extract time information
time_values <- time(b)

# Step 3: Turning strings into time
date_objects <- as.Date(time_values)

# Step 4: Get unique years
unique_years <- unique(format(date_objects, "%Y"))
#All years take too long, so we will only choose two years
unique_years <-c("2021", "2022")

# Step 5: Create an empty list to store mean stars for each year
mean_stars_list <- list()

Working with Temperature

This is where our programming skills come in handy from earlier lectures

# Step 6: Iterate over each unique year
for (year in unique_years) {
  cat("Calculating average for year", year, "\n")
  # Select dates for the current year
  selected_dates <- unique(date_objects[format(date_objects, "%Y") == year])
  # Find the index of the selected dates in the time_values
  selected_indices <- which(time_values %in% as.Date(selected_dates))
  # Extract raster values for the selected dates
  selected_rasters <- b[[selected_indices]]
  # Extract layer names
  layer_names <- names(selected_rasters)
  # Find indices of layers that start with "tmp"
  indices_tmp <- grep("^tmp", layer_names)
  # Extract temperature layers
  tmp_layers <- selected_rasters[[indices_tmp]]
  # Change raster names to temperature
  names(tmp_layers) <- selected_dates
  # Convert spatraster into raster
  stars_list <- lapply(tmp_layers, function(raster) st_as_stars(raster::raster(raster)))
  # Use lapply to apply st_apply to each stars object to calculate the mean
  mean_list <- lapply(stars_list, function(stars_obj) st_apply(stars_obj, c("x", "y"), mean))
  # Combine the list of results into a single stars object
  mean_stars <- do.call(stars::st_as_stars, mean_list)
  # Change the name of the raster
  names(mean_stars) <- paste("avg_tmp", year, sep = "_")
  # Add the mean stars to the list
  mean_stars_list[[year]] <- mean_stars
  cat("Average calculation for year", year, "completed\n\n")
}
Calculating average for year 2021 
Average calculation for year 2021 completed

Calculating average for year 2022 
Average calculation for year 2022 completed

Working with Temperature

This is where our programming skills come in handy from earlier lectures

# Combine all the mean stars into a single stars object
# Initialize an empty data frame
result_df <- data.frame()

# Loop through each year and calculate the mean
for (year in names(mean_stars_list)) {
  #Creating an object where we look through the years
  avg_tmp_col_name <- paste0("avg_tmp_", year)
  #Calculating mean and keeping track of the year
  mean_value <- mean(mean_stars_list[[year]][[avg_tmp_col_name]], na.rm = TRUE)
  #Saving everything into a dataframe
  result_df <- rbind(result_df, data.frame(year = as.numeric(year), average_temperature = mean_value))
}

# Writing the df to a dataframe
write.csv(result_df, "./data/output/average_temp_2021_2022.csv", row.names = FALSE)

Working with Temperature

So, we wrote our dataframe to a csv file.

Note that we did this dataframe for two years only: 2021 and 2022.

To complete the loop for all years, put a # in front of
unique_years <-c("2021", "2022")

This means that unique_years is an object made out of
unique(format(date_objects, "%Y"))

In other words, all years from 1901 until 2022. But this takes a long time (i.e. a few hours): multiply how long it took for 2021 and 2022 by 50.

To save you time, I have already done all the years on my machine. You can download the csv file with all the years.

Put it in the relevant folder

Working with Temperature

#write.csv(result_df, "./data/output/average_temp_1901_2022.csv", row.names = FALSE)
#Reading the file
average_temp<-read.csv("./data/output/average_temp_1901_2022.csv")
#Subsetting to below 2001
average_temp_1901_2000<-subset(average_temp, year<2001)
#Calculating the 1901-2000 average
average_temp_1901_2000<-mean(average_temp_1901_2000$average_temperature)
#Calculation the deviations from the 1901-2000 average
average_temp$dev_1901_2000_avg<-average_temp$average_temperature-average_temp_1901_2000

Working with Temperature

Creating a time trend

#Mapping the Temperature Deviaitonsover time
ggplot(average_temp, aes(x = year, y = dev_1901_2000_avg)) +
  geom_line() +
  geom_point() +
  theme_bw()+
  labs(title = "Time Trend of Average Temperature",
       x = "Year",
       y = "Average Temperature")

Working with Temperature

Creating a time trend

Working with Temperature

Creating a bar plot with the values

# Create a new column for bar color based on the condition
average_temp$bar_color <- ifelse(average_temp$dev_1901_2000_avg > 0, "red", "blue")

# Create the bar plot
ggplot(average_temp, aes(x = year, y = dev_1901_2000_avg, fill = bar_color)) +
  geom_bar(stat = "identity", color = "black") +
    scale_fill_manual(values = c("blue", "red"), guide="none") +

  theme_bw() +
  labs(title = "Time Trend of Average Temperature",
       x = "Year",
       y = "Difference from 1901-2000 avearge (degrees C)")

Working with Temperature

Creating a bar plot with the values

Working with Temperature

This is very close to the official figures produced by NOAA (National Ceneter for Environmental Information)

See: https://www.ncei.noaa.gov/access/monitoring/climate-at-a-glance/global/time-series/globe/land/1/11/1901-2022.

Raster Algebra

Without realizing, we have already covered an important component of working with rasters: raster algebra.

Raster algebra is the application of functions on a cell-by-cell basis, repeatedly for all pixels of one or more overlapping rasters.

Example:

Raster Algebra

Without realizing, we have already covered an important component of working with rasters: raster algebra.

Raster algebra is the application of functions on a cell-by-cell basis, repeatedly for all pixels of one or more overlapping rasters.

Example:

Raster Algebra

Without realizing, we have already covered an important component of working with rasters: raster algebra.

Raster algebra is the application of functions on a cell-by-cell basis, repeatedly for all pixels of one or more overlapping rasters.

Example:

Raster Algebra

Without realizing, we have already covered an important component of working with rasters: raster algebra.

Raster algebra is the application of functions on a cell-by-cell basis, repeatedly for all pixels of one or more overlapping rasters.

Example:

Raster Algebra

Without realizing, we have already covered an important component of working with rasters: raster algebra.

Raster algebra is the application of functions on a cell-by-cell basis, repeatedly for all pixels of one or more overlapping rasters.

Example:

Raster Algebra

Without realizing, we have already covered an important component of working with rasters: raster algebra.

Raster algebra is the application of functions on a cell-by-cell basis, repeatedly for all pixels of one or more overlapping rasters.

Example:

Raster Algebra

Raster Algebra Operations

  • Arithmetic: +, -, *, /
  • Logical: <, <=, >, >=, ==, !=, !

We can also apply the following operations to individual rasters: is.na, abs, round, sqrt, log, etc.

Raster Algebra Illustration

Let’s focus once again on Italy

#Raster Addition

We can add the two rasters together in the following way:

#Raster Average

But it probably makes sense to create an average of the two rasters

#Raster Average

As we showed previously, can make raster raplacements by using the dataframe operations

More Raster Algebra

  • As we saw previously, satellite luminosity is used as a proxy for economic activity

Luminosity in 2013

Working with Luminosity

These are the steps to have the data ready for plotting in ggplot

#Step1: Read country borders
borders <- st_read(dsn="./data/boundaries/all-country-boundaries.shp", quiet = TRUE)
#Step2: Creating a folder called
target_dir <- "./data/luminosity/F101992"
#Step3: Create the target directory if it doesn't exist
if (!dir.exists(target_dir)) {
  dir.create(target_dir)
}
#Step4: Unzip the file
untar("./data/luminosity/F101992.v4.tar", exdir = target_dir)
#Step5: Unzip the file further
R.utils::gunzip("./data/luminosity/F101992/F101992.v4b_web.avg_vis.tif.gz", remove = FALSE, overwrite=TRUE)
#Step6: Read the raster
f1992<-read_stars("./data/luminosity/F101992/F101992.v4b_web.avg_vis.tif")
#Step7: Make sure the two files have the same CRS
st_crs(f1992)<-st_crs(borders)
#Step8: Reduce the resolution of the raster
f1992b<-st_warp(f1992, st_as_stars(st_bbox(f1992), dx = 0.7))
#Step9: Rename the raster
names(f1992b)<-"lights_1992"

fig<-ggplot()+
  geom_stars(data=f1992b)+
  scale_fill_gradientn(colours=c("black","white"))+
  #scale_fill_viridis_c(option = "magma",begin = 0.1)+
  geom_sf(data=borders, linewidth = 0.1, fill = NA, color = "white", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

fig

Luminosity in 1992

Luminosity in 1992 in Europe

library(purrr)

#Identifying coordinates of Northern points
norway<-subset(borders, SOVEREIGNT == "Norway")
nc_geom <- st_geometry(norway)
list_points_norway<-nc_geom[[1]]
list_points_norway<-flatten(list_points_norway)
list_points_norway2<-as.data.frame(do.call(rbind, list_points_norway))
names(list_points_norway2)<-c("lon_x", "lat_y")
max_lat_y_eu<-max(list_points_norway2$lat_y)

#Identifying coordinates of Southern points
greece<-subset(borders, SOVEREIGNT == "Greece")
nc_geom <- st_geometry(greece)
list_points_greece<-nc_geom[[1]]
list_points_greece<-flatten(list_points_greece)
list_points_greece2<-as.data.frame(do.call(rbind, list_points_greece))
names(list_points_greece2)<-c("lon_x", "lat_y")
min_lat_y_eu<-min(list_points_greece2$lat_y)

#Identifying coordinates of Eastern points
ukraine<-subset(borders, SOVEREIGNT == "Ukraine")
nc_geom <- st_geometry(ukraine)
list_points_ukraine<-nc_geom[[1]]
list_points_ukraine<-flatten(list_points_ukraine)
list_points_ukraine2<-as.data.frame(do.call(rbind, list_points_ukraine))
names(list_points_ukraine2)<-c("lon_x", "lat_y")
max_lon_x_eu<-max(list_points_ukraine2$lon_x)

#Identifying coordinates of Western points
portugal<-subset(borders, SOVEREIGNT == "Portugal")
nc_geom <- st_geometry(portugal)
list_points_portugal<-nc_geom[[1]]
list_points_portugal<-flatten(list_points_portugal)
list_points_portugal2<-as.data.frame(do.call(rbind, list_points_portugal))
names(list_points_portugal2)<-c("lon_x", "lat_y")
min_lon_x_eu<-min(list_points_portugal2$lon_x)

Luminosity in 1992 in Europe

Mapping

fig<-ggplot()+
  geom_stars(data=f1992b)+
  scale_fill_gradientn(colours=c("black","white"))+
  geom_sf(data=borders, linewidth = 0.1, fill = NA, color = "white", alpha=0.5)+
  coord_sf(xlim = c(min_lon_x_eu-3, max_lon_x_eu+3), ylim = c(min_lat_y_eu-1, max_lat_y_eu-6), expand = FALSE)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fig

Luminosity in 1992 in Europe

Mapping

Luminosity in 1992 in the US

We can do the same about the US

min_lon_x_us<-(-130)
max_lon_x_us<-(-57)
min_lat_y_us<-(26)
max_lat_y_us<-(53)

fig<-ggplot()+
  geom_stars(data=f1992b)+
  scale_fill_gradientn(colours=c("black","white"))+
  geom_sf(data=borders, linewidth = 0.1, fill = NA, color = "white", alpha=0.5)+
  coord_sf(xlim = c(min_lon_x_us-3, max_lon_x_us+3), ylim = c(min_lat_y_us-3, max_lat_y_us+3), expand = FALSE)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fig

Luminosity in 1992 in the US

We can do the same about the US

Luminosity in 2013

Luminosity in 2013 in Europe

Luminosity in 2013 in the US

Comparing Luminosity

Comparing Luminosity

Positive Changes: 2013-1992>0

Comparing Luminosity

Here is how we do that:

#Step8: Reduce the resolution of the raster
f1992c<-st_warp(f1992, st_as_stars(st_bbox(f1992), dx = 0.5))
#Step9: Rename the raster
names(f1992c)<-"lights_1992"
#Step8: Reduce the resolution of the raster
f2013c<-st_warp(f2013, st_as_stars(st_bbox(f2013), dx = 0.5))
#Step9: Rename the raster
names(f2013c)<-"lights_2013"
#Step10: Turning to grid
f1992bpo<-st_as_sf(f1992c, as_points = FALSE, merge = FALSE)
f2013bpo<-st_as_sf(f2013c, as_points = FALSE, merge = FALSE)
#Step11: Spatially intersection 1992 and 2013
joined_data <- st_join(f1992bpo, f2013bpo, join = st_intersects)
#Switched off spherical geometry (s2)
sf::sf_use_s2(FALSE)
#Step12: Keeping only the ADMIN variable
borders2<-subset(borders, select = c("ADMIN"))
#Step13: Simplifying the world polygons
borders2<-st_simplify(borders2,  dTolerance = 0.005)
#Step14: Spatially intersecting borders and grids
joined_data_w<- st_join(joined_data, borders2, join = st_intersects)
#Step15: Removing seas
joined_data_w2<-subset(joined_data_w, !is.na(ADMIN))
#Step16: Calculating changes in luminosity
joined_data_w2$diff_lum<-joined_data_w2$lights_2013-joined_data_w2$lights_1992
joined_data2<-subset(joined_data_w2, diff_lum>0)
joined_data3<-subset(joined_data_w2, diff_lum<(-5))

Comparing Luminosity

Here is how we do that:

Positive Changes: 2013-1992>0

ggplot()+
  geom_sf(data=joined_data2, color=NA, aes(fill=diff_lum))+
  scale_fill_gradientn(colours=c("yellow", "red"), na.value = NA)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "black", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Comparing Luminosity

Here is how we do that:

Positive Changes: 2013-1992>0

Comparing Luminosity

Positive Changes in Europe: 2013-1992>0

ggplot()+
  geom_sf(data=joined_data2, color=NA, aes(fill=diff_lum))+
  scale_fill_gradientn(colours=c("yellow", "red"), na.value = NA)+
  geom_sf(data=borders, linewidth = 0.1, fill = NA, color = "black", alpha=0.5)+
    coord_sf(xlim = c(min_lon_x_eu-3, max_lon_x_eu+3), ylim = c(min_lat_y_eu-1, max_lat_y_eu-6), expand = FALSE)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Comparing Luminosity

Positive Changes in Europe: 2013-1992>0

Comparing Luminosity

Positive Changes in the US: 2013-1992>0

ggplot()+
  geom_sf(data=joined_data2, color=NA, aes(fill=diff_lum))+
  scale_fill_gradientn(colours=c("yellow", "red"), na.value = NA)+
  geom_sf(data=borders, linewidth = 0.1, fill = NA, color = "black", alpha=0.5)+
  coord_sf(xlim = c(min_lon_x_us-3, max_lon_x_us+3), ylim = c(min_lat_y_us-3, max_lat_y_us+3), expand = FALSE)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Comparing Luminosity

Positive Changes in the US: 2013-1992>0

Comparing Luminosity

Negative Changes: 2013-1992<0

ggplot()+
  geom_sf(data=joined_data3, color=NA, aes(fill=diff_lum))+
  scale_fill_gradientn(colours=c("yellow", "blue"), na.value = NA)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "black", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Comparing Luminosity

Negative Changes: 2013-1992<0

Comparing Luminosity

Negative Changes in Europe: 2013-1992<0

ggplot()+
  geom_sf(data=joined_data3, color=NA, aes(fill=diff_lum))+
  scale_fill_gradientn(colours=c("yellow", "blue"), na.value = NA)+
  geom_sf(data=borders, linewidth = 0.1, fill = NA, color = "black", alpha=0.5)+
    coord_sf(xlim = c(min_lon_x_eu-3, max_lon_x_eu+3), ylim = c(min_lat_y_eu-1, max_lat_y_eu-6), 
             expand = FALSE)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Comparing Luminosity

Negative Changes in Europe: 2013-1992<0

Comparing Luminosity

Negative Changes in the US: 2013-1992<0

ggplot()+
  geom_sf(data=joined_data3, color=NA, aes(fill=diff_lum))+
  scale_fill_gradientn(colours=c("yellow", "blue"), na.value = NA)+
  geom_sf(data=borders, linewidth = 0.1, fill = NA, color = "black", alpha=0.5)+
  coord_sf(xlim = c(min_lon_x_us-3, max_lon_x_us+3), ylim = c(min_lat_y_us-3, max_lat_y_us+3), expand = FALSE)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Comparing Luminosity

Negative Changes in the US: 2013-1992<0

Temperature

Temperature in 2019

library(terra)
library(sf)
#Step1: Reading the raster
b <- rast("./data/temp/cru_ts4.07.1901.2022.tmp.dat.nc")
#Step2: Turning it into stars object
stars_obj1 <- st_as_stars(raster::raster(b$tmp_1424))
#Step3: Renaming the raster
names(stars_obj1)<-time(b$tmp_1424)
#Step4: Ensuring the the same CRS
st_crs(stars_obj1)<-st_crs(border)

Temperature

Temperature in 2019

fig<-ggplot()+
  geom_stars(data=stars_obj1)+
   scale_fill_gradient2(
    low = 'blue',
    high = 'red',
    mid = 'yellow', 
    midpoint = 5,
    na.value = NA)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "white", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Temperature

Temperature in 2019

Temperature Difference

#Step1: Reduce the resolution of the raster
stars_1901<-st_warp(mean_stars_1901, st_as_stars(st_bbox(mean_stars_1901), dx = 1))
stars_2022<-st_warp(mean_stars_2022, st_as_stars(st_bbox(mean_stars_2022), dx = 1))

#Step2: Turning to grid
grid1901<-st_as_sf(stars_1901, as_points = FALSE, merge = FALSE)
grid2022<-st_as_sf(stars_2022, as_points = FALSE, merge = FALSE)

#Step3: Spatially intersection 1992 and 2013
joined_data <- st_join(grid1901, grid2022, join = st_intersects)

#Step4: Calculating the difference
joined_data$diff_temp<-joined_data$avg_tmp_2022-joined_data$avg_tmp_1901
joined_data2<-subset(joined_data, diff_temp>4)
joined_data3<-subset(joined_data, diff_temp<(-5))

Plotting Changes

Temperature difference between 2022 and 1901

fig<-ggplot()+
  geom_sf(data=joined_data, color=NA, aes(fill=diff_temp))+
  scale_fill_gradientn(colours=c("yellow", "blue"), na.value = NA)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "black", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Plotting Changes

Temperature difference between 2022 and 1901

Plotting Changes

Temperature difference between 2022 and 1901 > 4

fig<-ggplot()+
  geom_sf(data=joined_data2, color=NA, aes(fill=diff_temp))+
  scale_fill_gradientn(colours=c("yellow", "blue"), na.value = NA)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "black", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Plotting Changes

Temperature difference between 2022 and 1901 > 4

Popescu (JCU): Lecture 9